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In  this  paper  we  introduce  the  numerical  Laplace  transform,  a  local  time-frequency 
analysis  method  which  applies  to  causal  signals,  *The  numerical  Laplace  transform 
resolves  the  identity,  has  good  time-frequency  resolution,  and  adapts  resolution 
windows  according  to  the  time  delay.  The  numerical  Laplace  transform  is  equivalen : 
to  a  wavelet  transform  in  the  frequency  domain.  The  discretized  version  of  the 
numerical  Laplace  transform  is  invertible.  The  kernel  vectors  of  the  transform 
are  frame  vectors  that  are  nearly  tight  over  a  fairly  wide  range  of  parameters. 

We  demonstrate  this  with  several  numerical  experiments.  The  numerical  Laplace 
transform  resolves  a  causal  signal  onto  the  s-plane.  With  a  suitable  mapping, 
the  signal  is  resolved  into  the  frequency-damping  unit  disc. 
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Abstract 

In  this  paper  we  introduce  the  numerical  Laplace 
transform,  a  local  time- frequency  analysis  method 
which  applies  to  causal  signals.  The  numerical 
Laplace  transform  resolves  the  identity,  has  good  time- 
frequency  resolution,  and  adapts  resolution  windows 
according  to  the  time  delay.  The  numerical  Laplace 
transform  is  equivalent  to  a  wavelet  transform  in  the 
frequency  domain.  The  discretized  version  of  the  nu¬ 
merical  Laplace  transform  is  invertible.  The  kernel 
vectors  of  the  transform  are  frame  vectors  that  are 
nearly  tight  over  a  fairly  wide  range  of  parameters.  We 
demonstrate  this  with  several  numerical  experiments. 
The  numerical  Laplace  transform  resolves  a  causal  sig¬ 
nal  onto  s-plane.  Wiih  a  suitable  mapping,  the  signal 
is  resolved  into  the  frequency- damping  unit  disc. 


1  Introduction 

I  ho  literature  of  statistical  signal  processing  has 
been  dominated  over  the  past  twenty  years  by  work  on 
exact  and  subspace  methods  for  identifying  damped 
and  undamped  complex  exponential  modes  [4,  3]. 
Such  modes  arise  out  of  the  physics  of  lumped  linear 
systems,  so  they  naturally  carry  physical  information. 
Except  for  some  recent  extensions  [6],  the  topic  is  well 
developed. 

During  the  period  of  development  of  modal  anal¬ 
ysis,  there  has  been  renewed  interest  in  methods 
of  time-frequency  analysis  such  as  windowed  Fourier 
analysis  ,  \V  igner  distributions,  wavelet  decomposi¬ 
tions.  and  other  techniques  for  computing  infinite¬ 
dimensional  decompositions  of  signals.  In  a  sense, 
these  too  are  modal  decompositions,  for  they  decom¬ 
pose  a  signal  into  spectral  modes  that  are  localized  in 
time  and  frequency  or  time  and  scale.  However,  it  is 
usually  impossible  to  make  direct  links  between  the 


time-frequency  coefficients  and  physical  phenomena. 
Windowed  Fourier  analysis  does  not  help  because  the 
windowing  and  shifting  process  in  the  transform  mech¬ 
anism  simply  cut  the  signals  into  pieces.  It  is  therefore 
desired  to  develop  a  time-frequency  analysis  method 
which  bridges  the  gap  between  infinite- dimensional 
modal  decompositions  and  finite-dimensional  modal 
analysis. 


2  Continuous  Numerical  Laplace 
Transform 


We  start  from  a  wavelet  transform  in  the  frequency 
domain.  Define  the  mother  wavelet 


2” 


r(^) 


r(n)  (1  -|_  ji/)^ 


n  >  1 


(1) 


with  corresponding  time  domain  window  function 

(2)  - 

u{t)  is  the  unit  step  function,  n  G  N,  a  >  0 

Note  that  ^„(f)  has  unit  norm.  From  this  mother 
wavelet,  construct  the  basis  of  scaled  translates 


,  cr  >  0,  — oo  <  u;  <  oo 
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Then  define  the  complex  frequency  s  =  <r  +  jw  and 
write  the  scaled  translates  of  the  basis  as 


(3) 

(4) 

(5) 


s*  is  the  complex  conjugate  of  s.  The  inverse  Fourier 
transform  of  (jV;  s)  is 


3  Time-frequency  Resolution  of  the 
Continuous  Numerical  Transform 


n  —  l 

t  ^  € 


S*t 


(6) 


Define  a  wavelet  transform  of  an  analytic  signal  /(jV) 
as 

(W5;; /)(«)  =  (/(».  s))  (7) 

This  wavelet  transform  resolves  analytic  signals  onto 
the  s-plane.  The  wavelet  transform  resolves  the  iden¬ 
tity,  meaning  it  has  an  inverse.  Daubechies  [1]  calls 
it  an  isometry  from  the  Hardy  space  H^(  of  analytic 
signals)  onto  the  Bergman  space  of  analytic  functions 
on  the  upper  half  plane.  By  the  Parse val  identity, 
one  can  expect  a  time  domain  version  of  this  wavelet 
transform.  Now  we  are  ready  to  define  the  continuous 
numerical  Laplace  transform. 


The  continuous  numerical  Laplace  transform  has 
good  time-frequency  resolution.  One  way  to  see  this  is 
to  view  this  transform  as  a  variable  windowed  Fourier 
transform  [5]  with  2^„(crf)}  or  as 

the  family  of  window  functions.  Each  window  function 
has  centers  i*  and  lj*  and  radii  At  and  A^; 
in  the  time  domain  and  the  frequency  domain  respec¬ 
tively: 
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(13) 
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3 

(14) 

(15) 

A 

"  2Vn  -  2 

(16) 

Definition  1 

Let  f{i)  6  be  a  causal  signal.  Then  the 

nth  order  numerical  Laplace  transform  of  f{t)  written 

as  ('CtA„/)(s)  is 

(^tt„/)(s)  =  (/(<),  V'n(<;s)>  (8) 


where 

o'>0,u;GR,  n>l 


Therefore  the  resolution  cell  defined  by  the  window 
function  centered  at  (n/2(T,a;)  is 

2cr  2(7 

The  size  of  a  resolution  cell  is  thus 

(2A,)(2Aa,)  =  n>2 

y/n  — 2 

— ►  2,  as  n  — 00 


n 
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n  y/n 
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This  transform  resolves  causal  signals  onto  the  s- 
plane.  We  emphasize  that  this  is  not  quite  the  usual 
Laplace  transform,  because  the  kernel  is  essentially 
for  n  >  1,  not  Rewrite  the  kernel 

Vn(L  •*') 

ln{t:s)  =  (10) 

=  (11) 

One  sees  that  the  kernel  is  a  complex  expo¬ 

nential  with  "delay"  n  and  “damping”  (7. 

Tin*  inversf'  continuous  numerical  Laplace  trans¬ 
form  is  [5] 

Jo  J-CO 

(12) 


Note  that  this  size  varies  as  a  function  of  n.  When  ii 
increases  it  approaches  2  which  is  the  lower  bound  of 
the  Heisenberg  Uncertainty  Principle  which  can  only 
be  achieved  by  Gabor  transform  [2]. 

The  resolution  cells  are  not  rigid.  They  narrow  in 
time  resolution  when  damping  is  large  and  widen  when 
damping  is  small.  This  is  exactly  what  one  would 
hope  to  have  in  the  frequency-damping  modal  analy¬ 
sis.  Furthermore,  the  ratio  of  the  “center-time”  to  the 
“time  width”  is 

n/2(T  y/n 
2Vn/2(7  ^ 

which  is  indepent  of  the  damping  of  the  window  func¬ 
tion.  This  is  analogous  to  the  constant- Q  property 
in  wavelet  analysis.  A  plot  of  the  tiling  is  shown  in 
Figure  1. 


4  Discretized  Numerical  Laplace 
Transform  and  Frame  Bound  Esti¬ 
mation 


The  continuous  numerical  Laplace  transform  can  be 
discretized  to  a  lattice  in  the  s-plane,  namely  Smk  = 
-f  A  suitable  choice  is  am  =  a^  and  LJk  = 
a^kujQ  for  some  (Tq  >  0  and  a;o  >  0.  The  discretized 
numerical  Laplace  transform  is  defined  as  follows. 

Definition  2 

Let  f{t)  E.  be  a  causal  signal  and 

■>Pn{t;Smk)  -  Smk  =  <^0*  + 

ja^kiJo.ao  >  O.uq  >  0.  Then  the  nth  order  dis¬ 
cretized  numerical  Laplace  transform  of  f{t)  is 


where 


(£^J)im,k) 


(/(0>  ^m.fc)) 


(17) 


(18) 


m^k  E  Z 
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Tile  discretized  numerical  Laplace  transform  de¬ 
fined  in  the  Definition  2  is  invertible  only  when  the 
{(Pmk  =  V'n{t\Smk)]  Constitute  a  frame,  i.  e.  when 
there  exist  frame  bounds  0  <  A  <  B  such  that 

(19) 

mk 


for  all  /  in  the  underlying  Hilbert  space  [1].  Every 
function  /  in  the  underlying  Hilbert  space  can  then 
be  decomposed  and  reconstructed  as 

/  =  =  (20) 

mk  mk 


0)0 

A 

B 

B/A 

.25 

36.256 

36.262 

1.000 

.50 

18.127 

18.132 

1.000 

.75 

12.060 

12.112 

1.004 

1.00 

8.947 

9.182 

1.026 

1.25 

6.982 

7.522 

1.077 

1.50 

5.594 

6.492 

1.160 

(a) 


Wo 

A  ' 

B 

B/A 

.25 

13.019 

14.138 

1.083 

.50 

6.546 

7.092 

1.083 

.75 

4.364 

4.728 

1.083 

1.00 

3.223 

3.596 

1.116 

1.25 

2.001 

3.454 

1.726 

1.50 

0.325 

4.221 

12.986 

( 

Table  1:  Frame  bounds  for  (a)  numerical  Laplace 
frames  based  on  ips  and  (b)  wavelet  frames  based 
on  Maxican  hat  function  M{t)  =  2/\/37r“^/^(l  — 
^2)g-t  /2  dilation  prameter  =  2  in  both  cases. 


where  <j)mk  is  called  the  dual  frame  of  When 

^  ^  we  say  the  frame  is  tight  and  the  dual  frame 

is  just  the  frame  itself  scaled  by  the  inversion  of  the 
frame  bound.  In  this  case 

/  =  ^  ^  (A  <l>mk)  <t>mk  (21) 

mk 

When  a  frame  is  nearly  tight,  i.  e.  A  is  close  to  B, 
equation  (21)  serves  as  a  good  approximation. 

Since  «)}  define  a  continuous  wavelet  trans¬ 

form  ,  the  existence  of  numerical  Laplace  frames  is 
guaranteed  at  least  for  sufficiently  small  wo  [1].  More 
than  this,  it  turns  out  that  numerical  Laplace  frames 
are  nearly  tight  for  a  fairly  wide  range  of  choices  of  <to, 
cjo,  and  order  n.  Table  1  gives  the  frame  bound  esti¬ 
mations  for  numerical  Laplace  frames  of  order  3  and 
the  wavelet  frames  based  on  the  Maxican  hat  func¬ 
tion  [1].  (the  frame  bound  estimations  associated 
with  the  Maxican  hat  function  is  borrowed  from  [1]). 
Frame  bound  estimates  for  numerical  Laplace  frames 
for  other  choices  of  order  n  and  ao  are  also  given  in 
Table  2,  3,4,5. 


Wo 

A 

B 

B/A 

.25 

18.105 

18.154 

1.003 

.50 

9.050 

.75 

5.985 

6.102 

1.020 

1.00 

4.339 

1.25 

3.253 

■ilBBliM 

1.50 

2.473 

3.570 

1.443 

Table  2:  Frame  bounds  for  numerical  Laplace  frames 
based  on  ^5  with  (Tq  =  2 


Wo 

A 

B 

B/A 

.25 

12.004 

12.169 

1.014 

.50 

5.996 

6.090 

1.016 

.75 

3.917 

4.140 

1.057 

1.00 

2.745 

3.299 

1.202 

1.25 

1.956 

2.878 

1.471 

1.50 

1.394 

2.634 

1.889 

Table  3:  Frame  bounds  for  numerical  Laplace  frames 
based  on  V’7  with  <to  =  2 


Wo 

A 

B 

B/A 

.25 

22.678 

23.067 

1.018 

.50 

11.338 

11.539 

1.018 

.75 

7.540 

7.712 

1.023 

1.00 

5.578 

5.586 

1.051 

1.25 

4.321 

4.829 

1.118 

1.50 

3.425 

4.200 

1.226 

l'a[>|e  1:  Frame  bounds  for  numerical  Laplace  frames 
based  on  V3  with  (Tq  =  3 


Wo 

.4 

B 

B/A 

.25 

17.402 

18.855 

1.083 

.50 

8.700 

9.428 

1.084 

.75 

5.783 

6.302 

1.090 

1.00 

4.268 

4.796 

1.124 

1.25 

3.294 

3.957 

1.201 

1.50 

2.952 

3.451 

1.332 

Tablf  f)  Frame  bounds  for  numerical  Laplace  frames 
based  on  i  3  with  =  4 
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Figure  2:  Numerical  Laplace  frame  reconstruction  of 
7(t;3, 1)  with  n  =  5, AT  =  20, Ml  =  —5, M2  =  3,a;o  = 
0.75,(70  =  2, bound  zz  (A B)/2 

5  Numerical  Experiments 

Since  the  numerical  Laplace  frames  are  nearly  tight, 
the  tight  frame  decomposition  and  reconstruction  for¬ 
mula  (21)  can  be  used  as  a  good  approximation  with 
the  tight  frame  bound  approximated  by  {A  -f  B)I2 
of  the  numerical  Laplace  frames.  For  simplicity  we 
demonstrate  two  cases:  (1)  f{t)  =  7(^;3, 1)  for  a  nu¬ 
merical  Lapace  frame  of  order  5  and  (2)  f{t)  =  sin{t) 
for  a  numerical  Lapace  frame  of  order  7.  For  each 
case  we  present  the  experiment  results  with  two  plots. 
One  is  the  plot  of  the  original  signal  and  its  numer¬ 
ical  Laplace  reconstruction.  The  other  is  the  plot  of 
numerical  laplace  coefficients  in  gray  scales  over  the 
frequency-damping  plane,  while  the  former  demon¬ 
strates  the  good  approximation  of  tight  frame  expan¬ 
sions,  the  later  shows  how  the  energy  of  the  signal  is 
distributed  over  the  time-frequency  plane  in  numerical 
Laplace  frame  representations. 


6  Frequency-Damping  Resolution  of 
the  Unit  Disc 

Define  damping  coefficient  p  and  warped  frequency 
6  such  that 

P  =  e-\  pG(0,l)  (22) 

^  =  2tan”^a;,  ^G(“7r,  tt)  (23) 

Define  2  =  pe^^ .  Then,  via  this  mapping,  the  signal  is 
resolved  into  the  frequency-damping  unit  disc  in  the 
z-plane.  We  may  call  the  corresponding  transform  in 
z  variable  the  numerical  Z  transform  except  that  one 


NUfonn  witt)  f-gpdf(t;3,1).N-20.Ml-5.M2>3.wO-0.750000,»gO>2  000000 
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Figure  3:  The  numerical  Laplace  coefficients  in  gray 
scales  of  3, 1)  with  n  =  h^N  =  20, Ml  =  —5, M2  = 
3,a;o  =  0.75,0-0  =  2, bound  =  (vi  +  J5)/2 
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Figure  4;  Numerical  Laplace  frame  reconstruction  of 
7(^.3.  1)  with  n  =  7,N  =  20, Ml  =  -5,M2  =  4,a;o  = 
0. 75.0-0  =  2, bound  =  (>1  +  B)/2 


Figure  5:  The  numerical  Laplace  coefficients  in  gray 
scales  of  7(/; 3, 1)  with  n  =  7,A^  =  20, Ml  =  —5, M2  = 
4,a;o  =  0.75,0-0  =  2, bound  =  (^4  -f  B)/2 


should  be  reminded  that  this  transform  takes  contin¬ 
uous  signals  as  input. 

Figure  6  shows  the  lattice  after  the  mapping.  The 
locations  of  the  resolving  points  are  just  where  we 
would  expect  them  to  be:  finely  spaced  in  angular 
frequency  when  damping  is  small  and  widely  spaced 
when  damping  is  large.  This  should  serve  a  good 
model  for  frequency-damping  modal  analysis. 


Figure  6:  Resolving  points  of  the  numerical  Z  trans¬ 
form. 
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